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We study the variational solution of generic interacting fermionic lattice systems using fermionic 
Gaussian states and show that the process of gaussification, leading to a nonlinear closed equation 
of motion for the covariance matrix, is locally optimal in time by relating it to the time- dependent 
variational principle. By linearising our nonlinear equation of motion around the ground-state fixed 
point we describe a method to study low-lying excited states leading to a variational method to 
determine the dispersion relations of generic interacting fermionic lattice systems. This procedure 
is applied to study the attractive and repulsive Hubbard model on a two-dimensional lattice. 



The low-lying excited states of a quantum system are 
of central importance because they are, in contrast to 
the ground state, accessible in experiments. For example, 
spectral functions can be directly measured, revealing the 
particle content of a system. However, many theoretical 
models involve complicated many-body interactions, and 
an exact solution for their ground state and excitations 
are not known. In these cases, we must rely on numerical 
methods to find good approximations to the low-energy 
states of the system. 

In the case of one-dimensional lattice systems, the 
Density Matrix Renormalization group (DMRG) [U |2] 
has emerged as the most successful numerical method. 
By understanding the DMRG as an application of the 
variational method to the Matrix Product State (MPS) 
[3, 4 class, many generalizations have been recently pro- 
posed. This growing collection of powerful algorithms 
now allows for the determination of the ground state, 
low-lying eigenstates, and time-evolution of generic one- 
dimensional quantum spin systems O |6] . 

Comparatively less progress has been achieved in un- 
derstanding the low-energy physics of fermionic systems; 
fermions are the building blocks of all matter, and are 
thus central to many exciting effects in condensed matter 
physics, including superconductivity, superfluidity, and 
the Quantum Hall effect. While an MPS-based approach 
can be directly applied to one-dimensional fermionic sys- 
tems, via the Jordan Wigner transformation, this is no 
longer possible for higher-dimensional systems. Here, 
fermionic tensor networks, including fermionic Projected 
Entangled Pair States (PEPS) [7 and the Multiscale 
Entanglement Renormalisation Ansatz (MERA)[8 have 
been recently proposed for their study. However, so far, 
none of these approaches have been applied to describe 
the excitations of a fermionic quantum system. 

Another possibility to find ground-state approxima- 
tions of fermionic systems is to use other variational trial 
wave functions such as Slater determinants (Hartree Fock 
theory) or Gutzwiller projected wave functions. These 
have also been used to describe excitations in the past. 



in the framework of time-dependent Hartree Fock the- 
ory (see e.g. jHlllO]) and time-dependent Gutzwiller the- 
ory Recently, a novel technique based on the vari- 
ational class of fermionic Gaussian states (fGS)[12 has 
been proposed [13 . This method allows one to find ap- 
proximate solutions for ground- and thermal states, as 
well as the time evolution of interacting fermionic lattice 
systems in any dimension and geometry. An applica- 
tion to the 2d spinfull Fermi Hubbard mode has estab- 
lished that these algorithms are stable and efficient, and 
provided results in agreement with, and going beyond. 
Quantum Monte Carlo. The fGS class generalizes the 
basic building blocks of our understanding of fermionic 
matter, namely BCS-states for superfluid phases, and 
Slater determinants (Hartree Fock Theory) for Mott and 
fermionic spin states to the framework of generalized 
Hartree Fock Theory (gHFT)[14 , which is a strict su- 
perset of Hartree-Fock theory. As fGS can be described 
using a number of variational parameters scaling polyno- 
mially in the system size, it allows the efficient simulation 
of large systems. 

In this Letter, inspired by the utility of the fGS class 
for capturing a wide range of physically relevant phases, 
we take the next step and attempt to describe the low- 
lying excited states of interacting fermionic lattice sys- 
tems. Hence we begin by introducing the general model 
of interacting fermions we investigate in the framework of 
gHFT. Then we apply the Time Dependent Variational 
Principle (TDVP) [H |T5] to derive a locally (in time) 
optimal effective equation for the time evolution within 
the fGS variational class. An expansion of this equation 
around the ground-state solution thus leads to an ap- 
proximation of the low- lying excited states. We demon- 
strate the power of our approach by investigating the 
two-dimensional Hubbard model, both in the attractive 
and repulsive regimes, deriving the excitation spectrum, 
and discussing the nature of the excitations. 

We consider a system of M fermionic modes localised 
on a lattice in d dimensions. This is described by 
fermionic creation and annihilation operators obeying the 
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canonical anti-commutation relations (CAR), {a^,a^} = 
^k,h {o^kiCii} = 0- Such a system can be equivalently de- 
scribed by 2M real Majorana operators, Cj = <^ ] + % and 

Cj+M = - obeying the CAR {cj.Ck} = 28jk. 

Many systems of importance in condensed-matter physics 
are modeled by Hamiltonians involving two-body interac- 
tions. In the Majorana language, the most general such 
model is written as 



H = i^TklCkCi + ^ UklmnCkClCmCn 



(1) 



kl 



kin 



where T = -T^ G R^mm Ukimn G R is antisym- 
metric under the exchange of any two adjacent indices. 

Recently, in [13], the fCS variational class has been 
used to find approximate solutions for the ground- and 
thermal states of ([T]). The fCS class is defined as the set 
of all states which are exponentials of a quadratic form 
in the fermionic operators. Such states fulfill Wick's the- 
orem, and thus can be completely described on a single- 
particle level via the real and antisymmetric covariance 
matrix (CM) Vki = i/2([c/e, q]). For physical states the 
CM obeys the inequality < 1, while for pure states we 
have F^ = —1. Note that every pure fCS is the ground 
state of a quadratic Hamiltonian H = i^^ihkiCkCi^ 
where = — /i^ G R. Further, Gaussian states re- 
main Gaussian under time evolution according to any 
quadratic Hamiltonian, so that the time evolution can be 
formulated in terms of the CM alone as V{t) = 4[/i, F(t)]. 

We briefiy review how to approximate the time evo- 
lution of an interacting system of the form ([T]) using 
fCS. Since H generically includes nonquadratic interac- 
tions, any infinitesimal time step At takes us out of the 
fCS class. Thus, we have to project hack into the set 
of fCS after each time step. In p!3], this was done via 
the process of gaussification, where Wick's theorem is 
invoked after each time step to reexpress the system's 
state as an fCS. In this way, a closed evolution of the 
CM could be derived: t{t) = 4[h^^\r{t)),r{t)], where 
= T + 6tr2[/7F]. Thus time evolution can be 
formulated as the nonlinear evolution according to a 
quadratic but state- dependent Hamiltonian. In the fol- 
lowing, we use the Time- dependent Variational Principle 
to prove that gaussification, used ad hoc in [TT, is actu- 
ally locally optimal in time within the set of fGS. 

We briefiy recall the TDVP. Suppose we have some 
variational class of states {|?/^(x))|x G R^}. We 
aim to solve the time-dependent Schrodinger equation 
d/dt\ip{x.)) = —iH\ip{x.)) as best as possible while re- 
maining in our class. Unfortunately, in general, the vec- 
tor zi^|?/;(x)) is not an element of the tangent space, 
and we can only find an approximate solution. The 
optimal way to do this is to carry out the minimisa- 
tion miux \ \x^ dj\tlj{x.)) + zi^|?/^(x))||, were we have writ- 
ten l'^(x)) as a linear combination of vectors in the 
tangent space, and dj = d/dxj. If we perform this 




FIG. 1: Graphical representation of the TDVP within the 
variational manifold M of Gaussian states represented by 
state vectors |'0(x)). The time evolved vector iH\'ip{x.)) is 
in general not an element of the manifold and we have 
to find an optimal approximation within the variational man- 
ifold. Thus, we obtain the optimal path |'0(x)) within the 
manifold M. 



optimisation at each time step, and take the limit of 
infinitesimal step size, we arrive at the locally opti- 
mal evolution within our variational class. Note that 
this minimization is equivalent to solving for the Euler- 
Lagrange Equation 4t-^ — ^ = with Lagrangian 



£=(V(x)|(-i|-J?) |^(x)). 

We now show that gaussification is locally optimal by 
demonstrating its equivalence to the TDVP. To see this, 
note that any pure fermionic Gaussian state can be rep- 
resented as IV^g(^)) = Texp f Jq |0), 

where G = —G^ is a real matrix. Thus F and G 
are related via F = [G,F]. This relation is imposed 
in a Lagrangian using a Lagrange multiplier A, and 
we obtain for the TDVP Lagrangian the expression 
C = (V^(x)l {-^i-H) IV^(x)) = tr[F(G + h^'^T))] + 
tr[A(f - [G,F])], where h^^\r) = T + 3tr2[/7F]. The 
Euler-Lagrange equations of motion immediately lead to 
F = 4[/i^^\ F], which are exactly those arising from gaus- 
sification in Ref. [13]. Thus we obtain a locally optimal 
smooth path in our variational manifold, as depicted in 
Fig.[l] 

The optimality of the TDVP is now exploited to de- 
rive an approximation of the excitation spectrum within 
the fGS variational class. To this end we linearise 
our TDVP equation of motion around the variational 
ground state given by Fq, i.e., we write T{t) = Fq + 
eFi + O(e^), with real Fi = -Ff . [Note that Fq 
can be obtained via imaginary time evolution with re- 
spect to the quadratic, but state dependent, Hamilto- 
nian H = i^f. ih'^^^CkCi [Is]: starting from an arbi- 
trary pure state r(0) the evolution F(t) = O{t)r{0)O{t)^ 



where 0{t) = Texp 



£dt'2[r{t'),h^^\r{t'))] leads to 

Fq.] We can eliminate the time dependence by writing 
Ti{t) = e*'^^r(cj) and arrive at the following eigenvalue 
equation, which is the main result of this work: 



itotiito) = [/i(6)(ro),fi(w)] +6[tr[[/fi(u;)],ro]. (2) 
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While gives the general solution for the excitation 
spectrum within the fGS variational class, this matrix 
equation is, in general, hard to solve. The problem may 
be simplified by rewriting it as an eigenvalue equation 
via vectorization: 

iu\ri) = v\Ti), 

V = /i^^^ + 1 - 6(r(g)i + 1 (g)r)/7, (3) 

where V is a 4M^ x 4M^ matrix. Thus, the problem 
reduces to diagonalizing the matrix V, which can be nu- 
merically demanding for large systems. The problem is 
easier in the presence of additional symmetries, however: 
in the case of a translation invariant system the prob- 
lem simplifies considerably. In this case the matrix Fi 
is block-diagonal in Fourier space, i.e. Fi = 0j^Fi(k), 
with modes k and — k paired. Thus, for each momen- 
tum mode k, we only need to solve the reduced equation 
icj(k)|Fi(k)) = V(k)|Fi(k)). In this way we can extract 
the dispersion relation cj(k). Further, since the CM con- 
tains all information about a fGS, we can also infer the 
nature of the excitations by looking at Fi(k). 

For the remainder of this Letter we apply our ap- 
proach to the two-dimensional Hubbard model on a 
square lattice. The Hubbard model describes an interact- 
ing fermionic lattice system of particles with two internal 
spin states, a =t, |, 

^Hubb = ^ Q^x,a%,a + ^ ^xt^xj - M ^xg, 

(x,y),cr X x,cr 

(4) 

where x denotes a position on the lattice, (...) indicates 
a summation over nearest neighbors, and = <^]ccr^xcr 
is a particle number operator. We take the hopping pa- 
rameter t to be real and consider u positive (negative) in 
case of a repulsive (attractive) interaction. The chemical 
potential ja fixes the filling of the lattice. 

Despite its simple structure the Hubbard model allows 
for a wide range of physical phases, including Mott and 
spin ordered phases in the case of a repulsive interaction, 
and superfiuid phases in case of an attractive interaction. 
It is even believed that the doped Hubbard model at pos- 
itive u may provide a description of high-temperature su- 
perconductivity. However, unless we consider very spe- 
cial parameter regimes, an exact solution of the model 
is unknown, and despite an intense theoretical and nu- 
merical effort the precise structure of its phase diagram 
remains an open question. 

In the following we determine the excitation spectrum 
of the translation-invariant Hubbard model with periodic 
boundary conditions (PBC) in two dimensions within the 
fGS variational class using Eq. ([2|. To calculate V we 
need the CM of the ground state, Fq. For the transla- 
tion invariant case with PBC it has been shown in [M] 
that Fq can be obtained via a two-parameter optimiza- 



tion for aribitrary filling, both in the attractive and repul- 
sive regimes. We use this result to numerically determine 
Fq on a 31 X 31 lattice, transform into Fourier space, and 
are left, for each k, with the diagonalization of an 8 x 8 
matrix. 

Let's discuss our results for the Hubbard model in the 
attractive case for three exemplary sets of parameters 
{u, fi) = (-4, 1), (-2, 2), (-4, 3). In Table[l]we have sum- 
marised various properties of the ground states. We have 
calculated the filling n = N/ (2 • 31^) and the pairing per 

particle [16 , p = T.i,j,a,a' K4aSV')lV^. where N is the 
number of particles in the lattice. 





(-4,1) 


(-2,2) 


(-4,3) 


n 


0.83 


0.24 


0.17 


P 


0.028 


0.045 


0.13 



TABLE I: Filling n and pairing per particle p for the attractive 
Hubbard model on a 31 x 31 lattice. 

We consider configurations far from half-filling. This 
is because at half filling the Hubbard model has addi- 
tional symmetries leading to a highly degenerate ground 
state which severely complicates analysis of the excita- 
tions. The attractive Hubbard model supports superfiuid 
phases indicated by a non- vanishing pairing per particle. 
Thus, the ground state is gapped for all three sets of 
parameters. 

Next, we investigate the excitation spectra. The re- 
sults are depicted in Fig. [2] We see that for (li, ji) = 
(—4,1) and (li, /i) = (—4,3) the excitation spectrum is 
gapped, while it is gapless for {u, ja) = (—2,2). Even 
though the results shown in Fig. [2] seem to imply that 
there is only one possible excitation, a closer investigation 
reveals that for all three sets of parameters we have six 
non-trivial dispersion relations cjj(k), j = 1, . . . , 6, with 
ui (k) < a;2 (k) < — We find for all three sets of param- 
eters that UJ2 = 0J3 and uj^ = 6^5 = ujq. In Table [TT| we 
characterize these excitations via charge, spin and pairing 
order parameters C, S and A respectively. We find that 
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As 


s. 


St 
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(-4,1), (-4, 4) 


CJi, (JJ4-6 


<^2-3 










(-2,2) 




<^2-3 


C^4-5 


^4-5 






St = (4t^-H 


fe^), Afc = (4t4^)' = (4t«-fct) 
),Sz = (rifct - nki), C = (rifct + ^h) 



TABLE II: Classification of the excitations of the attractive 
Hubbard model. We distinguish charge (C), spin (S) and 
pairing (A) excitations and order uiCk) < Ct;2(k) < . . .. We 
find 6 non-trivial excitations in all three cases. Excitations 
ji < . . . < j2 with the same dispersion are labeled by ojj-^^-j^. 

For (u, /j) = (—2,2) we find that uoq has a spin and 
pairing order parameter of less than 1% compared to uj^-b 
and is thus set to zero in the table. 

for (li, /i) = (— 4, 1), (— 4, 3) the excitations are always 
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gapped, while they become gapless for (ix, /i) = (—2,2). 
Here, uq has a pairing Ak=o and a spin order Sz of less 
than 1% compared to the values for uj^-^. 
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FIG. 2: (Color online) Dispersion relation for the negative u 
Hubbard model 

Next, we discuss excitations in the repulsive Hub- 
bard model. Here, as predicted in the ground 
state is never paired, i.e. p = 0. Spin-ordered phases 
within a gHFT treatment are only predicted for half- 
filling, or large u. Thus, for our sets of parameters, 
(li, /i) = (4, 0), (4, —3), (4, —4) we only present the filling 
n in Table ini 





(4,0) 


(4, -3) 


(4,-4) 


n 


0.18 


0.70 


0.81 



TABLE ni: Filling n for the repulsive Hubbard model on a 
31 X 31 lattice. 



Now we discuss the excitation spectrum depicted in 
Fig. [3j We find ten non-trivial excitations for (li, /i) = 
(4, —3), (4, —4) [cjr-io have a fiat dispersion], and six 
for (ix,/i) = (4,0). Again, for each set of parameters 
there are excitations with different dispersion relations; 
we summarize the nature of the excitations in Table HVl 





Ak^o 


Ak 


As 


St 


(4, -3), (4, -4) 




C^3-4 
<^9-10 






(4,0) 


<^l-3, <^6 




(^1-3 


(^1-3 


As = (4^aLfc^), St = (al^a-ki) 



TABLE IV: Classification of the excitations of the attrac- 
tive Hubbard model. We distinguish spin (S) and pair (A) 

excitations and order (jUi(k) < (j02(^) < Excitations 

ji < . . . < j2 with the same dispersion are labeled by 
UJ7-1Q are the flat excitations near zero in Fig.|3] 

In contrast to the attractive Hubbard model, we find 
no charge excitations and only non-vanishing spin or- 
der St' The nature of the excitations for {u^ji) = 
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FIG. 3: (Color online) Dispersion relation for the positive u 
Hubbard model 



(4, —3), (4, —4) is identical, while we find a different struc- 
ture for {u, /J.) = (4,0). Here, the excitations with the 
highest energy is threefold degenerate, and we find no 
fiat excitations near zero energy. Note also, that for 
large negative chemical potential the excitations become 
gapped. 

In summary, we have presented how the the TDVP 
applied to the set of fermionic Gaussian states allows for 
a variational approach to determine dispersion relations 
and the nature of the excitations in interacting fermionic 
systems. For translational invariant systems, the aris- 
ing equations scale linearly in the system size, and can 
thus be applied to large systems in more than one dimen- 
sion, as we have demonstrated for the example of the 2d 
Hubbard model in both the attractive and the repulsive 
regimes on a 31 x 31 lattice. 
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